Apparatus and method for terrain model reproduction

ABSTRACT

The present invention relates to a method and an apparatus which reproduce a digital terrain model (DTM) from contour data with geomorphological consistency and natural features including fine folds. An initial DTM h o  is produced (step  101 ), and set h o  as an initial value (step  102 ). Then the operator T which smooths the elevational values along flowing water lines or the neighborhood of the lines is operated on the DTM h (step  103 ), further the operator B which sets the boundary values by substituting the contour data is operated on the DTM h (step  104 ). The number n of operating times of the steps  103  and  104  is checked (step  105 ). If n is less than preset number n max , return to the step  103  and if n reaches n max , go to the step  106  and the DTM h is output.

BACKGROUND OF THE INVENTION

1. Field of the Invention

The present invention relates to a method and apparatus which reproducea digital terrain model (DTM) with natural features including fine foldsfrom contour data.

2. Description of the Related Art

The simplest DTM from contour data is obtained by dividing a regionbetween contours and making the height of each part to be constant whichis the average of the elevational values of the enclosing contours. Itis called the stacking model. The information of the model is exactlyequivalent to that of the contour data. In the model every elevationalgap on contours and the flatness between contours become conspicuouswhere the density of contour is low. Hence it becomes a problem to varythe elevational data between contours and connect the contours by acontinuous surface. This is called the terrain reproduction problem fromcontour data. The necessary information to vary the elevational data isnot obtained directly from contour data. Therefore the problem is of anaddition of the information that some adequate information ought to beadded to contour data in order to obtain an appropriate height surfacedata. The problem can be divided into two partial problems; theframework of the addition and the content of the addition.

Most of conventional methods are classified into following four types.

(1) The profiles of a terrain are calculated for some directions usinginterpolation curves such as the spline, and averaged out with weights.

(2) Triangle patches are spanned throughout between contours, andelevational values are interpolated on the triangle patches.

(3) DTM is obtained by smoothing the initial model obtained by a simpleprocedure such as the stacking model by a two-dimensional low-passfilter.

(4) Regarding contours as a set of points, an estimating function isdeliberately chosen with a fitting surface. The surface is then fittedby minimizing the estimating function.

The outputs of above methods have some problems; (1) artificial stepsand ditches or starlike noises appear in DTM, (2) as shown in FIG. 5,the triangle patches 501 will remain as a peculiarity of the landform,(3) and (4) as shown in FIG. 6 the landform becomes an unusual roundedshape without minute folds and wrinkles.

Such problems occur in the conventional methods because physical andgeomorphological features of landforms are not given but only artificialand geometrical conditions as the content of the addition to contourdata. Also from the viewpoint of the framework of the addition, themethods are faced with some difficulties. The above (1), (2) and (3)methods are in principle of unified processes that the contents of theaddition are mostly determined according to the frameworks of theaddition. Thus it is almost impossible to add or cut a part of thecontents as the need arises.

In the method of (4), the content of the addition can be varied byexchanging the estimating function. However, the choice is small and itis very difficult to adopt a local requirement because the content isafter all given by a global-optimization problem.

SUMMARY OF THE INVENTION

A purpose of the present invention is to provide a method and anapparatus which reproduce a DTM from contour data with geomorphologicalconsistency and natural features including fine folds. It solves theproblems of geometrical noises and unnatural landform-features appearingin the result owing to merely geometrical and artificial techniques ofthe conventional methods.

To solve above problems in the present invention, the terrainreproduction problem is formulated as a boundary value problem for anoperator on a functional space defined on two-dimensional (2D) planewhere contours are set as the boundaries.

The content of the addition can be selected by the choice of theoperator and it is easy to adopt local requirements because theframework of the addition is the boundary value problem. Thus a widerange of physical-geomorphological features can be expressed.

Particularly, in an area where erosion by rain dominates the landformation, it is most natural to give the geomorphological feature thatthe elevational values are smoothed along flowing water lines to thecontent of the addition. It is fundamental in the terrain reproductionproblem on such an area. Then in the present invention, an operatorwhich smooths the elevational values along flowing water lines orneighborhoods of the lines is introduced and it makes possible to giveat least the feature to a DTM.

Some other features of rain erosion may be also given to a DTM byconstructing an operator to which corresponding operations are added.For example, an effect of accumulation can be adopted into a DTM with anoperator which behaves like the Laplace operator in flat areas.

In this application, an operator which smooths the elevational valuesalong flowing water lines or the neighborhood of the lines means anoperator at least including this smoothing operation.

As a result, a method and an apparatus of the present invention solvethe problems of geometrical noises and unnatural landform-featuresresulting from the conventional methods and it makes possible toreproduce a DTM from contour data with geomorphological consistency andnatural features including fine folds.

BRIEF DESCRIPTION OF THE DRAWINGS

FIG. 1 is a block diagram showing an apparatus for executing the DTMreproduction in Embodiment 1 or Embodiment 2 of the present invention.

FIG. 2 is a flowchart of Embodiments of the present invention.

FIG. 3 is a schematic view explaining an operator in Embodiment 1 of thepresent invention.

FIG. 4 is a schematic profile showing the discrepancy between the twointerpolation methods in Embodiment 1 of the present invention.

FIG. 5 is a schematic view showing the remaining triangle texture of thelandform by a conventional method using triangle patches.

FIG. 6 is a schematic view showing the unusual rounded shape of thelandform by a conventional method using a two-dimensional low-passfilter or a fitting surface by minimization.

FIG. 7 is a schematic view showing the visualization of a DTM bythree-dimensional projection in Embodiment 4 of the present invention.

DESCRIPTION OF THE PREFERRED EMBODIMENTS

The present invention will be explained below by embodiments.

First, setting a notation, the framework of the addition in theembodiments of the present invention will be expressed.

An arrow “A←B” means the procedure that a content of variable B issubstituted into variable A. Let a terrain reproducing domain be D. Forthe simplification, only a square domain is considered. The domain D isquantized by dividing D into N×N square blocks of same size and eachblock is expressed by a pair of integers (x,y) (x,y=0,1,2, . . . ,N−1).It should be noticed that the above quantization is introduced only forsimplifying the explanation. Other quantizations such as triangles withunequal size can be used.

The contour lines are one-dimensional subsets of D. They are expressedby C₁, C₂, . . . , and their elevational values by v₁, v₂, . . .respectively.

A DTM is expressed by a function of two integer variables (x,y). A valueof the function at (x,y) is the elevational value of a correspondinglattice point. Namely, it is expressed by a 2D array h(x,y), forexample.

An operator on a functional space defined on a 2D plane is atransformation of a 2D array into a 2D array. The operation of anoperator T transforming a 2D array h into a 2D array h′ is denoted byh′=T(h). The composite transformation of two operators S and T,h′=S(T(h)) is denoted by S·T.

An example of boundary value problem for an operator is that an equationof h including an operator, h=T(h) is solved for h under the boundarycondition that h(x,y)=v for every point (x,y) of C which is aone-dimensional subset of D.

A boundary value problem for an operator is to solve h=T(h), forinstance, which is an equation of h with an operator T for h under theboundary condition that h(x,y)=v for every point (x,y) of C which is anone-dimensional subset of D.

Under the above notation, one of the frameworks of the addition in thepresent invention is formulated by the following boundary value problem.

Let an operator which smooths the elevational values along flowing waterlines or the neighborhood of the lines be T. As a terrain model issmoothed, it finally arrives at the equilibrium state where no furthersmoothing proceeds.

Then the DTM h is expressed by the solution of the boundary valueproblem that it satisfies the equation h=T(h) under the condition thath(x,y)=v_(i) at every point (x,y) of C_(i) (i=1, 2, . . . )respectively.

The equation of this type is adopted because the solution is easilyobtained. The method of the present invention can be executed by someother formulation such as T′(h)=0, derived from the differentconsideration.

If an adequate initial value h_(o) is given and the operator T is stablein the neighborhood of the initial value, the boundary value problem canbe solved by iterative method using the operator B which sets theboundaries by substituting the contour data.

The operation B(h)=h′ is defined by the following rule oftransformation. If a point (x,y) is included in C_(i) (i=1,2, . . . ),then h′(x,y)=v_(i), otherwise h′(x,y)=h(x,y).

The operator S is defined by S=B·T. It is known that if h is calculatedby h=S·S· . . . ·S(h_(o)), operating S on h_(o) iteratively, h becomesan approximate solution for the boundary value problem.

Even when the iterative solution does not converge, an approximatesolution can be obtained asymptotically.

Since the equation itself is only an phenomenological model, it ismeaningless to require much of the strictness of the solution.

An iterative method other than the above can also be used.

There is the case that shapes of adjacent contours are noncorrelativefor some reasons such as two contours apart too much from each other.Even such the case, it is sometimes possible to obtain a natural resultby proceeding a modificatory filtering that is an iterative operation ofonly T (excluding B) to the obtained approximate solution h for severaltimes. This is a method that gives the content of the addition possessedby the operator priority over the information of contour data (It is aninvention corresponding to claim 6).

The hardware costitution of an apparatus of an embodiment will beexplained referring to the block diagram in FIG. 1.

The apparatus 1 of the embodiment has a memory device 2, a processor 3,an input device 4, an output device 5 and a control device 6, and theyare connected mutually by bus lines 7 and 8.

The memory device 2 is provided with:

a region for storing contour data, i.e., the contour data store unit;

a region for storing a program to produce the initial DTM, i.e., theinitial DTM production program unit;

a region for storing a program to operate the operator which smoothselevational values of the input DTM along flowing water lines or theneighborhood of the lines, i.e., the smoothing program unit;

a region for storing a program to set the boundary conditions, i.e., theboundary condition setting program unit;

a region for storing the DTM generated by the apparatus, i.e., the DTMstore unit; and

a region for storing a control program such as OS, i.e., the controlprogram unit.

The processor 3 is a CPU and other related wares. The input device 4constitutes of for example a digitizer, a mouse, a keyboard, a numericalfile, and possibly a light pen. The output device 5 is for example anumerical file or a relay of network. The control device 6 controls therespective devices for executing the program.

Next, the elementary action of the apparatus of the embodiment isexplained with reference to the flowchart of FIG. 2.

First, the initial DTM h_(o) is produced (step 101), and set the h_(o)as an initial value (step 102).

Then the operator T which smooths the elevational values along flowingwater lines or the neighborhood of the lines is operated on the DTM h(step 103), further the operator B which sets the boundary values bysubstituting the contour data is operated on the DTM h (step 104). Thenumber n of operating times of the steps 103 and 104 is checked (step105). If n is less than preset number n_(max), return to the step 103and if n reaches n_(max), go to the step 106 and the DTM h is output.

In the flowchart of FIG. 2, the first means which produces the initialDTM from contour data consists of the step 101 and the step 102, and thesecond means which smooths the elevational values along flowing waterlines or the neighborhood of the lines consists of the steps from 103 to106.

By selections of the operator and extensions of the constitution of theapparatus including the two means, various kinds of embodiments of thepresent invention will be possible.

Since the first means is common in every embodiment, it is explainedhere. The definitions of operators and the extensions of constitution ofthe apparatuses will be explained in the following embodiments.

As the first means, for example, the succeeding two procedures can beadopted. The stacking model is produced from contour data and it is setas the initial terrain model of the iterative solution method.

The second procedure is trivial, thus, the first is shown. A stackingmodel is easily made by an area-painting technique of computer graphics.One point is picked up in an area enclosed by contours and the point isset as the seed of the area-painting. Then the area is painted with thecolor expressing the elevational value, which can be determined from thecontours.

Other than the stacking model, some terrain models without conspicuousnoises, for example, an interpolated model with triangular patches canbe used for the initial terrain model.

EMBODIMENT 1

In a first embodiment, all the calculations at each point including thedetermination of flowing water lines and the smoothing of elevationvalues are executed on the point and its eight-points neighbourhood.Since the calculations are local, the integrator part of the flowingwater lines is dispensable and the implementation becomes easy.

The operation of the smoothing operator T on DTM h at a point (x,y) isdefined by the followings. The explanation is made referring to FIG. 3.

(1) From the elevational values of a point (x,y) and its eight-pointsneighborhood, the gradient vector grad h(x,y) is determined. Allelevational values of the points in the neighborhood should be used tocalculate the gradient vector.

(2) In the case that the gradient vector is not null on the point (x,y)(such a point is shown by P₁ in FIG. 3). The smoothing operator T_(o) istaken action. The operation h′=T_(o)(h) is defined as follows:

1) The gradient vector (indicated by 301 in FIG. 3) is quantized intothe closest one of four directions which are horizontal, vertical, 45degree upper-right and 45 degree upper-left. The gradient vector 301 inFIG. 3 is quantized as the 45 degree upper-right direction.

2) The point P₁(x,y) and two points P₁, P₂ from its eight-pointsneighbourhood (U₁, in FIG. 3) on which the line of the quantizeddirection passes are picked up and a set of these three points isconsidered as the local flowing water line 302 in FIG. 3.

3) Some weights, which total sum is one, is given to execute thesmoothing. The elevational values on the local flowing water line areaveraged with the weights, and this gives the new elevational valueh′(x,y) on the point P₁:(x,y). A set of weights whose elements are 0 forP₁ and ½ for P₂ and P₃ in FIG. 3 is given as one of the examples.

(3) In the case that the gradient vector is null (such a point is shownby P₄ in FIG. 3).

To extend a flowing water line (c in FIG. 3), the elevational values onthe point P₄ and its eight-points neighborhood (U₂ in FIG. 3) areaveraged with weights whose total sum is 1. The weights are set toexecute the two-dimensional smoothing. This procedure gives the newelevational value. A set of weights whose elements are ¼ for the fourshaded points around P₄ in FIG. 3 and 0 for the other points is given asone of the examples.

In the case of (2), the smoothing means is not only restricted to theabove but may be a filter directly depending on the gradient vector or anon-linear filter.

In the case of (3), the smoothing means can also be one of non-linearfilters. For example, a non-linear operator defined by

h′(x,y)=(h(x−1,y)+h(x+1,y))·|h _(x)|/4(|h _(x) |+|h_(y)|)+(h(x,y−1)+h(x,y+1))·|h _(y)|/4(|h _(x) |+|h _(y)|),

where h_(x)=h(x+1,y)−h(x−1,y) and h_(y)=h(x,y+1)−h(x,y−1), can be usedin (3). If a smooth DTM is given as the initial value, the operator thatonly consists of (1) and (2) may be used because the above case (3)rarely occurs for such the initial value. If the operator T isimplemented in the flowchart shown in FIG. 2 and the operations areexecuted, then the DTM is obtained. The present embodiment solves thedifficulty that geometrical noises and unnatural features appear in theoutput.

EMBODIMENT 2

In a second embodiment, the calculational area of the operator isextended to execute the smoothing on somewhat global area, and anintegrator becomes necessary to obtain the flowing water line withenough length. Further, the degree of the smoothing is decreased toobtain clear ridges and valleys.

The operation of the smoothing operator T on DTM h is defined below.

(1) The gradient vector field, grad h of the DTM h, is determined.

(2) From the vector field, discrete integral calculation is executed toobtain the flowing water line. The line consists of consecutive latticepoints that passes the point (x,y). Its parametrization is defined bythe following.

c(t)=(c ₁(t), c ₂(t)) (c(0)=(x,y), c(t)≠c(t+1))

It should satisfy the condition that a inner product (c(t)−c(t−1),c(t+1)−c(t)) is positive.

(3) In the case that the flowing water line does not terminate in anarea necessary for the smoothing. The smoothing operator T₁ whoseoperation h′=T₁(h) is defined by

h′(x,y)=Σw(c(t)−(x,y))h(c(t))/Σw(c(t)−(x,y))

is operated. Here, Σ expresses the sum for every t which satisfies thecondition: w(c(t)−(x,y))≠0. For example, the weight w is given by

w(x,y)=exp (−(x ² +y ²)/(2σ²)) (for −3σ<x,y<3σ, σ: constant),

w(x,y)=0 (otherwise).

The example is Gaussian distribution type. One of the simplest operatorsis shown here, and any other operator such as a non-linear filter can beused.

(4) In the case other than (3)

1) It is decided whether a point (x,y) belongs to a ridge-valley area ornot.

Here, a ridge-valley area denotes an area where a ridge or a valleyexists. The decision, for example, is executed by the followingprocedure: Variables a₊, a⁻, b₊ and b⁻ is introduced as

a ₊ =h(x+1,y)−h(x,y), a ⁻ =h(x−1,y)−h(x,y),

b ₊ =h(x,y+1)−h(x,y), b⁻ =h(x,y−1)−h(x,y).

A point is in a ridge-valley area if a₊·a⁻>0 or b₊·b⁻>0 satisfies andotherwise not. The procedure is simply one of the examples and a moreglobal procedure can be used for the decision.

2) In the case that a point (x,y) belongs to the ridge-valley area.

The smoothing operator T₂ whose operation h′=T₂(h) is defined by

h′(x,y)=(h(x,y)+w _(d)(x,y)·h _(L)(x,y))/(1+w _(d)(x,y))

is operated. Here, h_(L)(x,y) is defined as

h _(L)(x,y)=(h(x−1, y)+h(x+1,y)+h(x,y−1)+h(x,y+1))/4

and the w_(d)(x,y) is a weight depending on the degree of ridge orvalley on a point (x,y). For example, the weight is given by a functionof the Laplacian h (physically speaking, the amount of water flow into apoint). Various other techniques can be considered to define T₂. Forexample, W_(d) can be set as an increasing function which depends on thedistance from a point where the gradient vector vanishes and some 2Dnon-linear filters can be also used instead of the h_(L).

3) In the case other than 2).

The smoothing operator whose operation is defined byh′(x,y)=h_(L)(x,y)(h_(L) is defined above) is operated.

In the above (4), a simpler constitution can be considered whichexecutes the smoothing procedure of the above 3) without theridge-valley area decision. In this case, it is regarded as aconstitution simply replacing T_(o) of Embodiment 1 with T₁.

The above operator T is implemented in the flowchart in FIG. 2 and theoperations are executed to obtain the DTM. The present embodiment solvesthe problems in an area with low-dense contours, keeping the ridges andvalleys from being flattened and sustaining the sharpness of terrainsurface. As the result, it can reproduce the DTM from the contour datawith geomorphological consistency and natural features including finefolds.

EMBODIMENT 3

In a third embodiment, an extended constitution for the fast calculationwill be shown. It consists of an interpolation method and the aboveiterative operations with several unit sizes of quantization.

When the iterative method is used in solving the present boundary valueproblem, the number of the iteration necessary to reach the equilibriumstate of the smoothing is approximately proportional to the horizontalinterval of contours in general.

For example, if the unit size of the quantization is changed four timesas large as the original unit size and each side of the domain D isdivided into N/4, the number of blocks in each interval of contoursbecomes ¼ of the original number. The number of the iteration is thusreduced to ¼ of the original number. The number of lattice points whichare used in the calculation is also reduced to ¼×¼, namely {fraction(1/16)}. Hence, the amount of the calculation required to obtain theresult is reduced to {fraction (1/64)} comparing with that of theoriginal. The unit size of the quantization is now changed two times aslarge as the original unit size and a DTM is calculated by magnifyingthe obtained result through interpolation on a lattice whose size isN/2×N/2. This gives a much better approximation of the equilibrium statethan a stacking model. Therefore, the number of the iterations is verylimited at the step of the quantization. A DTM in the quantization ofthe original size is obtained by recursively executing the similarprocedure which magnifies the result through interpolation.

An example of the constitution of the present embodiment is as follows.

(1) The scale of the input contour data is reduced to ¼ in both sides,and a ¼ scale DTM is made by the apparatus of Embodiment 2. (It shouldbe noticed that every lattice in this constitution expresses theoriginal domain D with different quantizations.)

(2) The above result is magnified two times through interpolation on alattice with N/2×N/2 points.

(3) The above magnified result is regarded as the initial model and thescale of the input contour data is reduced to ½ in both sides. A ½ scaleDTM is made by the second means of the apparatus of Embodiment 2, usingthese as the boundary values.

(4) The above result is magnified two times through interpolation on theoriginal lattice.

(5) The above magnified result is regarded as the initial model. Theinput contour data is used as the boundary values and an original scaleDTM is made by the second means of the apparatus of Embodiment 2.

In the above example, the scale ratio of 2 is only used either for thereduction or the magnification at each step but other ratio can be used.Even it can be changed at every step of the quantizations. It is alsonot necessary to take ¼ as the starting scale ratio.

It is not suitable to calculate an elevational value of the interpolatedpoint from that of the nearest adjacent points by using linearinterpolation such as a plane fitting. The ridges and the valleys areneglected and the flatness is enlarged by such the interpolation. It isdifficult to restore the lost unevenness from the second means and theridges and the valleys are being crushed.

The unevenness of terrain can, however, be interpolated by a cubicfitting with such the least-square method or Lagrange interpolation. Aschematic profile of terrain is shown in FIG. 4 to explain theinterpolation, where the axis of abscissa represents the x-axis and theaxis of ordinate represents that of elevational values. The points P₁,P₂, P₃ and P₄ have known elevational values and P is a point whereelevational value is interpolated. The solid line c is the profile of alinearly interpolated plane and the broken line c′ is the profile of aninterpolated curved surface in higher order. The elevational value atthe point P is h for the linear interpolation and h′ for the higherorder interpolation. The discrepancy occurs between them as shown in thefigure. The elevational value h′ gives a summit whereas h simply aplateau. It is therefore better to use a higher order interpolation forthe magnification of DTM through interpolation.

The above method is known as a fast solution for partial differentialequations. It has been applied to such a boundary problem as Laplaceequation and has an effect several times faster than the ordinary methodon speed of calculation. Nevertheless, its effectiveness has not beenfully utilized in the conventional use.

It is quite clear from its constitution that it is most effective inobtaining a rough solution when the magnitude of frequency of thesolution gradually decreases as the frequency increases. In theconventional problem, the solution does not have such structure socalled fractal (or self-similar) structure, and further the strictnessof result is required for most of time.

On the contrary, these factors become advantages for a terrainreproduction problem. In recent years, it is elucidated that a terrainmodel has fractal structure from the success of fractal science ingeomorphology, and the strictness of result is not substantiallyrequired for the reproduction problem.

According to the test using contour data of several actual terrains, itis possible to accelerate the calculation speed by several decadal timesto about hundred times. Thus if it is compared with that of theconventional ordinary problems, the effect is more than ten times. Theremarkable effect must be realized that it is more than a known methodsimply applied to the new problem.

By using this fast calculation method, the number of iterations of thesmoothing operator can be tremendously reduced and this suppresses thesmoothing on areas of ridge or valley. Accordingly, a DTM is obtainedwith sharp and clear ridges and valleys.

EMBODIMENT 4

In a fourth embodiment, a display output process which visualize a DTMis explained referring to FIG. 7. The process of three-dimensionalprojection is here given as an example. It comprises in the followingsteps:

(1) A DTM is produced by the apparatus of Embodiment 3. The modificatoryfilter in claim 6, which uses the smoothing operator T of Embodiment 2,is applied to the DTM.

(2) An imaginary observation point (P_(o)), an imaginary light source(L) and an imaginary projection plane (H) are set first. Then each pointof the DTM (P_(D)) is projected to the projection plane. At the sametime, a tone of each point of the DTM on the projection plane (P_(H)) isdetermined by at least one from among the following elements;

a) the position of the point from the observation point,

b) the elevational value of the point,

c) the gradient of the surface at the point,

d) the ray direction at the point, and

e) the irradiance at the point.

As the projection method, a known method such as z-buffer or ray-tracingcan be used.

Further, a DTM can be also visualized by directly determining a tone ofeach point of a DTM by at least one from among b) to e) without animaginary observation point and an imaginary projection plane.

(3) The obtained toned figure is shown on a screen or printed out.

According to Embodiment 4, fine folds, ridges and valleys clearlyappear, and the natural terrain is successfully visualized as a whole.

While the preferred embodiments of the invention has been described,such description is for illustrative purpose only, and it is to beunderstood that changes and variations may be made without departingfrom the spirit or scope of the following claims.

What is claimed is:
 1. A terrain model production apparatus comprising:a first means for setting an initial terrain model calculated as afunction defined on two-dimensional plane from contour data; and asecond means for producing a terrain model by solving a boundary valueproblem including an operator which smoothes elevational values of aterrain model along flowing water lines or their neighborhoods with saidcontour data as boundary values by an iterative method with said initialterrain model.
 2. A terrain model production apparatus according toclaim 1, wherein said operator includes an integration means for atleast locally calculating the flowing water lines by integrating thegradient vector field of a terrain model, which is calculated beforehandand a smoothing means for smoothing unevenness of the elevational valueson each flowing water line of said terrain model along the respectiveflowing water line or its neighborhoods.
 3. A terrain model productionapparatus according to claim 1, wherein said operator includes adetection means for finding a ridge area and a valley area of a terrainmodel and a smoothing means for operating two-dimensionally in theneighborhoods of a point Where said gradient vector becomes zero, saidsmoothing means reducing the degree of the smoothing in an area found bysaid detection means.
 4. A terrain model production apparatus accordingto claim 1 further comprising: a third means for producing a terrainmodel of low precision from said first means and second means by settinga quantization unit in a defined domain of the terrain model larger thanrequired precision; a fourth means for producing the terrain model byapplying said second means to the initial value set by changing saidquantization unit smaller and correspondingly magnifying said terrainmodel through interpolation; and a fifth means for recursively producingthe terrain model by iteratively using the fourth means till reachingthe required precision.
 5. An apparatus according to claim 4, whereinsaid magnification through interpolation is a higher order interpolationusing the elevational values of larger area than the nearestneighbourhood of the interpolation points.
 6. A terrain modelmodification apparatus comprising: means for storing a terrain model asan initial value; and a modificatory filter operable once or recursivelyon the stored terrain model which smoothes elevational values alongflowing water lines or their neighborhoods.
 7. A terrain modelproduction method comprising: a first step of setting an initial terrainmodel calculated as a function defined on two-dimensional plane fromcontour data; and a second step of producing a terrain model by solvinga boundary value problem including an operator which smootheselevational values of a terrain model along flowing water lines or theirneighborhoods with said contour data as boundary values by an iterativemethod with said initial terrain model.
 8. A terrain model visualizationmethod comprising: a first step of producing a terrain model by themethod of claim 7; a second step of determining a tone of a point of theterrain model by setting at least one imaginary light source and usingat least one element from among elements which are the elevational valueof the point, the gradient vector of the terrain surface at the point,the ray direction at the point and the irradiance at the point; and athird step of visualizing a tone determined by the second step.
 9. Arecording medium in which a visualized data produced by a method ofclaim 8 is stored.
 10. A method for printing a terrain model, comprisingthe steps of: setting a toned figure produced by a method of claim 8;and printing the toned figure.
 11. A terrain model print produced by amethod of claim
 10. 12. A projected terrain model visualization methodcomprising: a first step of producing a terrain model by the method ofclaim 7; a second step of projecting a point of the terrain model withsetting an imaginary observation point, an imaginary projection planeand at least one imaginary light source; a third step of determining atone of a point of the terrain model by using at least one element fromamong elements which are the elevational value of the point, theposition of the point from the observation point, the gradient vecto ofthe terrain surface at the point, the ray direction at the point and theirradiation at the point; and a fourth step of visualizing a tonedetermined by the third step.
 13. A recording medium in which avisualized data produced by a method of claim 10 is stored.
 14. A methodfor printing a terrain model, comprising the steps of: setting a tonedfigure produced by a method of claim 12; and printing the toned figure.15. A terrain model print produced by a method of claim
 14. 16. Arecording medium in which a terrain model produced by a method of claim7 is stored.
 17. A terrain model production method comprising: a firststep of setting contour data as boundary values; and a second step ofcalculating a terrain model as a function defined on two-dimensionalplane by at least simulating rain erosion under a boundary condition setby the first step.
 18. A recording medium in which a terrain modelproduced by a method of claim 17 is stored.
 19. In a terrain modelingmethod operable on a computer for reproducing a terrain model withgeomorphological consistency and natural features including folds, theimprovement comprising: a first step of setting a terrain model as aninitial value; and a second step of operating an operator once orrecursively which smoothes elevational values of the terrain model alongflowing water lines or their neighborhoods.
 20. A recording medium inwhich a computer program is stored, said program being provided forexecuting: calculating a terrain model as a function defined on atwo-dimensional plane from contour data and setting said model as aninitial terrain model; and producing a terrain model by solving aboundary value problem including an operator which smoothes elevationalvalues of a terrain model along following water lines or theirneighborhoods with contours as the boundary values by an iterativemethod with said initial terrain model.